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ABSTRACT 

In the framework of the Lindblad theory for open quantum systems the damping of 
the harmonic oscillator is studied. A generalization of the fundamental constraints on 
quantum mechanical diffusion coefficients which appear in the master equation for the 
damped quantum oscillator is presented; the Schrodinger and Heisenberg representations 
of the Lindblad equation are given explicitly. On the basis of these representations it 
is shown that various master equations for the damped quantum oscillator used in the 
literature are particular cases of the Lindblad equation and that the majority of these 
equations are not satisfying the constraints on quantum mechanical diffusion coefficients. 
Analytical expressions for the first two moments of coordinate and momentum are also 
obtained by using the characteristic function of the Lindblad master equation. The master 
equation is transformed into Fokker-Planck equations for quasiprobability distributions. A 
comparative study is made for the Glauber P representation, the antinormal ordering Q 
representation and the Wigner W representation. It is proven that the variances for the 
damped harmonic oscillator found with these representations are the same. By solving the 
Fokker-Planck equations in the steady state, it is shown that the quasiprobability distribu- 
tions are two-dimensional Gaussians with widths determined by the diffusion coefficients. 
The density matrix is represented via a generating function, which is obtained by solving 
a time-dependent linear partial differential equation derived from the master equation. 
Illustrative examples for specific initial conditions of the density matrix are provided. 

1. Introduction 

In the last two decades, more and more interest arose about the problem of dissipation 
in quantum mechanics, i.e. the consistent description of open quantum systems [1-4]. The 
quantum description of dissipation is important in many different areas of physics. In 
quantum optics we mention the quantum theory of lasers and photon detection. There 
are some directions in the theory of atomic nucleus in which dissipative processes play 
a basic role. In this sense we mention the nuclear fission, giant resonances and deep 
inelastic collisions of heavy ions. Dissipative processes often occur also in many body or 
field-theoretical systems. 

The irreversible, dissipative behaviour of the vast majority of physical phenomena 
comes into an evident contradiction with the reversible nature of our basic models. The 
very restrictive principles of conservative and isolated systems are unable to deal with more 
complicated situations which are determined by the features of open systems. 

The fundamental quantum dynamical laws are of the reversible type. The dynamics of 
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a closed system is governed by the Hamiltonian that represents its total energy and which 
is a constant of motion. In this way the paradox of irreversibility arises: the reversibility 
of microscopic dynamics contrasting with the irreversibility of the macroscopic behaviour 
we are trying to deduce from it. 

One way to solve this paradox of irreversibility is to use models to which Hamiltonian 
dynamics and Liouville's theorem do not apply, but the irreversible behaviour is clearly 
present even in the microscopic dynamical description. The reason for replacing Hamilto- 
nian dynamics and Liouville's theorem is that no system is truly isolated, being subject 
to uncontrollable random influences from outside. For this reason these models are called 
open systems. There are two ways of treating quantitatively their interaction with the 
outside. One is to introduce specific stochastic assumptions to simulate this interaction, 
the other is to treat them according to the usual laws of dynamics, by regarding the open 
system as a subsystem of a larger system which is closed (i.e. which obeys the usual laws 
of dynamics, with a well-defined Hamiltonian). The dissipation arises in general from the 
subsystem interactions with this larger system, often reffered to as the reservoir or bath. 
The first of these two approaches has been used for the study of steady-state transport 
processes, in systems obeying classical mechanics. The second of the two approaches has 
been mainly used in quantum mechanics. The main general result [1, 5-7] is that under 
certain conditions the time evolution of an open system can be described by a dynamical 
semigroup ^t{t > 0). For a closed finite system with Liouville operator the evolution 
operator is not restricted to nonnegative t. The importance of the dynamical semigroup 
concept is that it generalizes the evolution operator to open systems, for which there is no 
proper Liouville operator and no $t for negative t. The mathematical theory of dynamical 
semigroups has been developed in [1, 8-12]. 

The quantum mechanics of the unidimensional damped harmonic oscillator represents 
a fundamental theoretical problem with applications in different domains of quantum op- 
tics, solid state physics, molecular and nuclear physics. In the present paper the quantum 
harmonic oscillator is treated in the Lindblad axiomatic formalism of quantum dynamical 
semigroups. 

In Sect. 2 the notion of the quantum dynamical semigroup is defined using the concept 
of a completely positive map [10]. The Lindblad formalism replaces the dynamical group 
uniquely determined by its generator, which is the Hamiltonian operator of the system, by 
the completely positive dynamical semigroup with bounded generators. Then the general 
form of Markovian quantum mechanical master equation is given. 

In Sect. 3 we give the fundamental constraints on quantum mechanical diffusion co- 
efficients which appear in the corresponding master equations [17]. The Schrodinger and 
Heisenberg representations of the Lindblad equation are given explicitly. On the basis of 
these representations it is shown that various master equations for the damped quantum 
oscillator used in the literature for the description of the damped collective modes in deep 
inelastic collisions or in quantum optics are particular cases of the Lindblad equation and 
that the majority of these equations are not satisfying the constraints on quantum me- 
chanical diffusion coefficients. Explicit expressions of the mean values and variances are 
also given [17,18]. 

In Sect. 4 we solve the master equation with the characteristic function [19]. This 
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function is found as a solution of a corresponding partial differential equation. By this 
method one can derive explicit formulae for the centroids and variances and, in general, 
for moments of any order. 

In Sect. 5 we explore the applicability of quasiprobability distributions to the Lind- 
blad theory [22] . The methods of quasiprobabilities have provided technical tools of great 
power for the statistical description of microscopic systems formulated in terms of the 
density operator [58]. The first quasiprobability distribution was the one introduced by 
Wigner [43] in a quantum-mechanical context. In quantum optics the P representation, 
introduced by Glauber [44,45,50], provided many practical applications of quasiprobabil- 
ities. The development of quantum-mechanical master equations was combined with the 
Glauber P representation to give a Fokker-Planck equation for the laser [47,48]. The 
master equation of the one-dimensional damped harmonic oscillator is transformed into 
Fokker-Planck equations for the Glauber P, antinormal Q and Wigner W quasiprobability 
distributions associated with the density operator. The resulting equations are solved by 
standard methods and observables directly calculated as correlations of these distribution 
functions. We solve also the Fokker-Planck equations for the steady state and show that 
variances found from the P, Q and W distributions are the same [22] . 

In Sect. 6 we study the time evolution of the density matrix that follows from the 
master equation of the damped harmonic oscillator [20,21]. We calculate the physically 
relevant solutions of the master equation by applying the method of generating function. 
This means that we represent the density matrix with a generating function which is the 
solution of a time-dependent partial differential equation of second order, derived from the 
master equation of the damped harmonic oscillator. We discuss stationary solutions of 
the generating function and derive the Bose-Einstein density matrix as example. Then, 
formulas for the time development of the density matrix are presented and illustrative 
examples for specific initial conditions provided. The same method of generating function 
was already used by Jang [39] who studied the damping of a collective degree of freedom 
coupled to a Bosonic reservoir at finite temperature with a second order RPA master 
equation in the collective subspace. 

The conclusions are given in Sect. 7. 

2. Lindblad theory for open quantum systems 

The standard quantum mechanics is Hamiltonian. The time evolution of a closed 
physical system is given by a dynamical group Ut which is uniquely determined by its 
generator H, which is the Hamiltonian operator of the system. The action of the dynamical 
group Ut on any density matrix p from the set D(H) of all density matrices of the quantum 
system, whose corresponding Hilbert space is H, is defined by 

p(t) = C/,(p) = e-i^Ve^^* 

for all t G (— cxD, cxd). We remind that, according to von Neumann, density operators 
p e D(H) are trace class {Trp < oo), self-adjoint (p^ = p), positive (p > 0) operators 
with Trp = 1. All these properties are conserved by the time evolution defined by Ut. 
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In the case of open quantum systems the main dificulty consists of finding such time 
evolutions $t for density operators p{t) = which preserve these von Neumann con- 

ditions for all times. From this requirement it follows that $t must have the following 
properties: 

(i) $t(Aipi + A2P2) = Ai$t(pi) + A2$t(p2); Ai, A2 > with Ai + A2 = 1, 
{ii) $t(p+) = Mp)+, 
{in) $t(p) > 0, 
(iv) Tr$t(p) = 1. 

But these conditions are not restrictive enough in order to give a complete description 
of the mappings $f as in the case of the time evolutions Ut for closed systems. Even 
in the last case one has to impose other restrictions to Uf, namely, it must be a group 
Ut+s = UfUs- Also, it is evident that in this case Uo{p) = p and Ut{p) ^ p in the trace 
norm when t ^ 0. For the dual group Ut acting on the observables A e B(H), i.e. on the 
bounded operators on H, 



Also, Ut{A) — > A ultraweakely when t — > and Ut is an ultraweakely continuous mapping 
[1,7,9,10,13]. These mappings have a strong positivity property called complete positivity: 



Because the detailed physically plausible conditions on the systems, which correspond 
to these properties are not known, it is much more convenient to adopt an axiomatic point 
of view which is based mainly on the simplicity and the succes of physical applications. 
Accordingly [1,7,9,10,13] it is convenient to suppose that the time evolutions ^t for open 
systems are not very different from the time evolutions for closed systems. The simplest 
dynamics $t which introduces a preferred direction in time, which is characteristic for 
dissipative processes, is that in which the group condition is replaced by the semigroup 
condition [6,7,11,12,14] 




J2B+Ut{A+Aj)Bj > 0,A,B, e B(H). 



The duality condition 




(2.1) 



and 



(2.2) 



are equivalent. Also the conditions 



(2.3) 
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ultraweakely when t — > and 

Mp) P 

in the trace norm when t — > 0, are equivalent. For the semigroups with the properties 
(2.2), (2.3) and 

A > ^ ^t{A) > 0, 

it is well known that there exists a (generally bounded) mapping L-the generator of 
$t is uniquely determined by L. The dual generator of the dual semigroup is denoted 
by L: 

Tr{L{p)A) = Tr{pL{A)). 
The evolution equations by which L{L) determine uniquely ^t{^t) are 



dt 

and 



= (2.4) 



^ = £(i.(^)), (2.5) 

respectively, in the Schrodinger and Heisenberg picture. These equations replace in the 
case of open systems the von-Neumann-Liouville equations 

and _ 
respectively. 

For any applications eqs.(2.4) and (2.5) are only useful if the detailed structure of the 
generator L{L) is known and can be related to the concrete properties of the open systems, 
which are described by such equations. 

Such a structural theorem was obtained by Lindblad [10] for the class of dynamical 
semigroups $t_which are completely positive and norm continuous. For such semigroups 
the generator L is bounded. In many applications the generator is unbounded. 

A bounded mapping, L : B(H) B(H) which satisfies L{I) = 0, L{A+) = L{A)+ 

and _ _ _ 

L{A+A) - L{A+)A - A+L{A) > 

is called dissipative. The 2-positivity property of the completely positive mapping 

^t{A+A) > $,(A+)$,(A), 

with equality at t = 0, implies that L is dissipative. Lindblad [10] has shown that con- 
versely, the dissipativity of L implies that is 2-positive. L is called completely dissipative 
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if all trivial extensions of L are dissipative. Lindblad has also shown that there exists a 

one-to-one correspondence between the completely positive norm continuous semigroups 
$t and completely dissipative generators L. The structural theorem of Lindblad gives the 
most general form of a completely dissipative mapping L [10,11]: 

Theorem: L is completely dissipative and ultraweakely continuous if and only if it is 
of the form 

LiA) = i[H,A] + ^ + [^^' ^]^^-), (2.6) 

j 

where V^,T,jV;Vj e B(H), H e B(H),.^.. 

The dual generator on the state space (Schrodinger picture) is of the form 

L{P) = ~[H,P] + ^J2(\-^jP^^j'^ + [V„pV+]). (2.7) 

j 

Eqs.(2.4) and (2.7) give an explicit form for the most general time-homogeneous quantum 
mechanical Markovian master equation with a bounded Liouville operator. 

Talkner [16] has shown that the assumption of a semigroup dynamics is only applicable 
in the limit of weak coupling of the subsystem with its environment, i.e. for long relaxation 
times. 

We should like to mention that all Markovian master equations found in the literature 
are of this form after some rearrangement of terms, even for unbounded generators. 

It is also an empirical fact that for many physically interesting situations the time 
evolutions $t drive the system toward a unique final state p(oo) = limt_>oo ^t(p(0)) for all 
p(0) e D(H). 

The evolution equations of Lindblad, being operator equations, the problem of finding 
their solutions is, in general, rather difficult. In cases when the equations are exactly 
solvable, these solutions give complete informations about the studied problem - they 
permit the calculation of expectation values of the observables at any moment. 



3. Master equations for damped quantum harmonic oscillator 

In this Section the case of damped quantum harmonic oscillator is considered in the 
spirit of the ideas presented in the previous Section. The basic assumption is that the 
general form (2.7) of a bounded mapping L given by Lindblad theorem [10] is also valid 
for an unbounded completely dissipative mapping L: 

Hp) = ~[H,p] + + [Vj,pv;]). (3.1) 

j 

This assumption gives one of the simplest way to construct an appropriate model for this 
quantum dissipative system. Another simple condition imposed to the operators Vj^ 
is that they are functions of the basic observables of the one- dimensional quantum me- 
chanical system q and p with [q, p] = ihl, where / is the identity operator on H of such 
kind that the obtained model is exactly solvable. A precise version for this last condition is 
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that linear spaces spanned by the first degree (respectively second degree) noncommutative 
polynomials in p and q are invariant to the action of the completely dissipative mapping 
L. This condition implies [15] that Vj are at most the first degree polynomials in p and q 
and H is at most a second degree polynomial in p and q. 

Beacause in the linear space of the first degree polynomials in p and q the operators 
p and q give a basis, there exist only two C-linear independent operators Vi, V2 which can 
be written in the form 

Vi = ttip + biq, i = 1,2 

with ai,bi = 1,2 complex numbers [15]. The constant term is omitted because its contri- 
bution to the generator L is equivalent to terms in H linear in p and q which for simplicity 
are chosen to be zero. Then H is chosen of the form 

H = Ho + ^{pq + qp), Ho = ^p' + ^q'. (3.2) 

With these choices the Markovian master equation can be written: 

dp % % % 

= -j\Ho,p\ - + fj.)[q, pp + pp] + ^i>^ - p)[p, pq + qp]- 

-^[q. [q.p]] - ^\p, \p,p]] + ^{[q, \p,p]] + [p, hp]])- (3.3) 

n n n 

Here we used the notations: 

i=l,2 3 = 1,2 j=l,2 i=l,2 

where Dpp,Dqq and Dpq are the diffusion coefficients and A the friction constant. They 
satisfy the following fundamental constraints as shown in [17]: 

i) Dpp > 

ii) Dqq > (3.4) 

in) DppDqq - Dpq'' > X'h'/4. 

Introducing the annihilation and creation operators 

-.{^/mu}q + ^^p), 



1 i 

a"*" = —={^mwq -. p), (3-5) 



obeying the commutation relation [a, a"*"] = 1, we have 

HQ^hu;{a+a+^) (3.6) 
and the master equation has the form 

-£ = -{Di - ii){pa+a+ - a+pa+) + ^{Di + n){a+a+p - a+pa+) + 



+ -(D2 — A — iuj){a^ pa — paa^) + -(-D2 + A + iu)){apa^ — a^ap) + /i.e., (3-7) 

where 

ft ma; 

L»2 = ^(ma;i^gg + ^). (3.8) 

In the hterature, equations of this kind are encountered in concrete theoretical models 
for the description of different physical phenomena in quantum optics, the damping of 
collective modes in deep inelastic collisions of heavy ions or in the quantum mechanical 
description of the dissipation for the one-dimensional harmonic oscillator. In the following 
we show that all these master equations are particular cases of the Lindblad equation 
and that the majority of these equations are not satisfying the constraints on quantum 
mechanical diffusion coefficients, and therefore the uncertainty principle is violated. 

1) The Dekker master equation for the damped quantum harmonic oscillator [4,23-26] 
supplemented with the fundamental constraints (3.4) obtained in [23] from the condition 
that the time evolution of this master equation does not violate the uncertainty principle 
at any time, is a particular case of the Lindblad master equation (3.7) when p, = \. 

2) The quantum master equation considered in [27,28] by Hofmann et al. for treating 
the charge equilibration process as a collective high frequency mode is a particular case of 
the Lindblad master equation (3.3) if A = 7(c<;)/2m = Dqq — 0, Dpp — 7(a;)T*(u;), Dpq = 
0, but the fundamental constraints (3.4) are not satisfied. 

3) For the quantum master equation considered in [29] for the description of heavy ion 
collisions we have X = n = 7/2, Dpp = D, Dqq = 0, Dpq = Dgp = —d/2 and consequently 
the fundamental constraints are not fulfilled. 

4) In [30], Spina and Weidenmiiller considered two kinds of master equations I and 
II for describing the damping of collective modes in deep inelastic collisions of heavy 
ions. Eq.I can be obtained from eq.(3.3) by replacing Hq with Hq — ^Amujq^ + f{t)Q ^-^d 
setting A = /U = T/2,Dpp = D/2, Dqq = and Dpq = Dqp = B/2. Then the constraints 
(3.4) are not satisfied. Eq.II is obtained from (3.3) by putting Hq — {l/2)A^^mujq^ — 
(l/2ma;)A^V + fit)q and = F^^ = A,// = Q.Dpp = D'//2,Dqq = Dj{/2,Dpq = 
and the last condition (3.4) is satisfied for all values of the parameters. 

5) The master equation for the density operator of the electromagnetic field mode 
coupled to a squeezed bath [31,32] can be obtained from the master equation (3.7) if we 
set 

/. = 0,A = 7,^(A-^^-#^) = -iV,^(^^-#^ + 2z%) = M. 
2A n nmuj 2A n nmuj n 

6) The master equation for the density operator of a harmonic oscillator coupled to an 
environment of harmonic oscillators considered in [33-36] is a particular case of the master 
equation (3.7) if we put 

A = ^ = -,,o„ = o,c„ = oi(A-^) = -n 
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and the fundamental constraints (3.4) are not fulfilled. 

7) The master equation written in [37] for different models of correlated- emission 
lasers can also be obtained from the master equation (3.7) by putting 

111 1 

-{Di + //) = A4, -{Di - //) = A3, -{D2 + X + iu;) = A2, -{D2-X- iu;) = Ai. 

8) Two master equations were introduced by Jang in [38,39], where the nuclear dis- 
sipative pocess is described as the damping of a collective degree of freedom coupled to 
a bosonic reservoir at finite temperature. The resulting RPA master equation within the 
observed collective subspace is derived in a purely dynamical way. The master equation 
written in [38] in the resonant approximation (rotating-wave approximation) can be ob- 
tained as a particular case of the Lindblad master equation (3.7). For this one has to 
set 

= m^uj^Dgg, Dpg = ii = 0, ^ = (2 < n > +l)r, A = -, 

where < n > is the average number of the RPA collective phonons at thermal equilibrium 
and F is the width (friction parameter). The fundamental constraints (3.4) are fulfilled in 
this case. 

The master equation derived recently [39] in order to extend the calculations carried 
out in [33] with the before-mentioned master equation in the resonant approximation, can 
also be obtained as a particular case of the master equation (3.7) by taking 

Dqq = Dpq = 0, Dpp = ^1^(2 < n > +l)r, = A = ^ 

or D2 = Di = {2 < n > -|-l)F/2 and in this case the fundamental constraints (3.4) are not 
fulfilled. 

The following notations will be used: 

a,{t)^Tr{p{t)q), 

ap{t) = Tr{p{t)p), 
aqq = Tr{p{t)q^)-al{t), (3.9) 
(Tpp = Tr{p{t)p^) - cTp(t), 

apq{t) = Tr{p{tf-^^)-ap{t)aq{t). 
In the Heisenberg picture the master equation has the following symmetric form: 

=L($,(A)) = '-[HoMA)] - 4(A + /.)([$t(A),g]p + #,(A),g]) + 
- + [$t(A),p]g) - ^[g, [g, $,(A)]]- 
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^'^''■[p, [p, MA)]] + %([P, h MA)]] + [q, \p, MA)]]). 



Denoting by A any selfadjoint operator we have 

aAit) = Tr{p{t)A),aAAit) = Tr{p{t)A^) - al{t). 

It follows that 

= TrL{p{t))A = Trp{t)L{A) (3.10) 

and 

= TrL{p{t))A^ - 2^^aA{t) = Trp{t)L{A') - 2aA{t)Trp{t)L{A). (3.11) 

An important consequence of the precise version of solvability condition formulated 
at the beginning of the present Section is the fact that when A is put equal to p or g in 
(3.10) and (3.11), then dap{t)/dt and daq{t)/dt are functions only of (7p{t) and (7q{t) and 
dapp{t) /dt, daqq{t) /dt and dapq{t)/dt are functions only of app{t),aqq{t) and apq{t). This 
fact allows an immediate determination of the functions of time crp(t), crq(t), app{t), aqq{t), 
(Tpq{t). The results are the following: 

^ = -(A-^)a,(t) + la,(t), 

= -mu^aq{t) - (A + p)ap{t) (3.12) 

= -2(A - p)aqq{t) + ^^pS) + '^Dqq, 

^ = -2(A + p)(Tpp{t) - 2muj'^apq{t) + 2Dpp, (3.13) 

= -mUj'^aqqit) + — CTpp(t) - 2Xapq{t) + 2Dpg. 

All equations considered in various papers in connection with damping of collective modes 
in deep inelastic collisions are obtained as particular cases of Eqs. (3.13), as we already 
mentioned before. 

The integration of Eqs. (3. 12) is straightforward. There are two cases: a) p, > oj 
(overdamped) and b) p < u> (underdamped). If S{t) denotes the vector 

<^<iit) 
ap{t) 

and M the 2x2 matrix 

— {\ — p) 1/m 



and 



M 



—muj —{X-\-p) 
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then (3.12) becomes 

^^ = MS{t). (3.14) 

tit' 

Now M can be written as M = N~^FN with F a diagonal matrix. It foUows that the 
solution of (3.14) is 

S{t) = N-^e^'NS{0). 
In the case a) with the notation v"^ — jj? — uP' the matrices N ^ and F are given by 

N 



-{X + u) 
-(X-u) 

■ sinh zvt cosh vt — ^ sinh z^t 7 ' 
aq{t) = e~'^*{(coshi't+ - sinhi/t)crq(0) + — smhi'tap{0)), 

2 

ap{t) = e~^*(-^^ sinhi/ta-q(0) + (coshz^t - ^ sinhzyt)a-p(O)). (3.15) 

If A > z^, then a"q(oo) = (Jp(cxo) = 0. If A < z^, then a"g(oo) = a"p(oo) = oo. In the case b) 
with the notation il^ = o;^ — /x^, the matrices N, and F are given by- 



Then 



I.e., 



N 



Then 



I.e., 



/_(A + zf2) 
' -{X-iQ))- 



1 X , / cos + T^Ot ^oSinOt 



sin cos Vtt — sin fit / ' 



u b If L 

2 

ap{t) = e-^\-^^^smntaq{0) + (cos lit - ^ sinl]t)c7p(0)) (3.16) 
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and CTq (oo) = (Tp(oo) = 0. 

In order to integrate Eqs.(3.13), it is convenient to consider the vector 

(^pq{t) J 

Then the system of equations (3.13) can be written in the form 

dX{t) 



dt 

where R is the foUowing 3x3 matrix 



RX{t) + D, 



-2{\-ii) 2uj 

R=\ -2(A + //) -2a; 

—uj uj — 2A 

and D is the foUowing vector 

(2mujDqq 
— Dr,r> 
2Dpg 

Then there exists a matrix T with property = I where / is the identity matrix and a 
diagonal matrix K such that R = TKT. From this it follows that 

X{t) = {Te^*T)X{0) + T(e^* - I)K-^TD. (3.17) 

In the overdamped case (/x > ui) the matrices T and K are given by 

^ / IJ, + v jj, — v 2lo 
T = — \ — V 11 + u 2u> 
\ —uj —ijj —2\i 

and 

(-2{\-v) 

K=\ -2{\ + v) 

\ -2A 

with = \^ — uj^ . 

In the underdamped case (// < a;) the matrices T and K are given by 

^ / jj, + iQ jji — iVt 2u3 
jjL — iVL iJL-\- iVL 2uj 



2iO. 



\ —UJ —U! —2fi 
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and 

/-2(A-zO) 

K ^ i -2{X + in) 

\ -2A 

with = u;^ — //^. 

From (3.17) it follows that 

X(oo) = -(TK-^T)D = -R-^D (3.18) 

(in the overdamped case the restriction A > is necessary). Then Eq.(3.17) can be written 
in the form 

X{t) = {Te^*T){X{0) - X{oo)) + X{oo). (3.19) 



Also 



and 



dt 

dX{t) 
dt 



= {TKe'^*T)(X(0) - X(oo)) = R(X{t) - X(oo)) 



t=o = {TKT){X{0) - X{oo)) = R{X{0) - X{oo)). 



The formula (3.18) is remarcable because it gives a very simple connection between the 
asymptotic values of (7qq{t) , (7pp{t) , crpq{t) and the diffusion coefficients Dqq, Dpp, Dpg. As an 
immediate consequence of (3.18) this connection is the same for both cases, underdamped 
and overdamped, and has the following explicit form: 

= 2(ma;)W + a;^-^^) ^^"^^^'^'^^^ + + 
+iij'^Dpp + 2rrujj'^{\ + ii)Dpq), 
(Tpp{oo) = 2A(A2 +L2 _ ^2) {{mwfuj^Dqq+{2X{X-fi)+uj^)Dpp-2mw\X-fi)Dpq), (3.20) 

(7pq(oo) = 2y^A(A^ +co^- fi'^) ^ p){mwfDqq + (A - ii)Dpp + 2m(A^ - ^i^)Dpq). 

These relations show that the asymptotic values aqq{po)^(7pp{oo)^apq{oo) do not depend 
on the initial values crgg(O), crpp(O), apq{Q). In other words, 

_^ / 2A(A + //) + a;^ 2uj{X + y) 

^~' = AU\2^ 2 2S\ 2A(A-//)+a;2 _2a;(A - /i) 

4A(A2 + a;2-/.2) 1 _ 2(A2-/,2) 



Conversely, if the relations D = — i?X(oo) are considered, i.e., 

/ 2mujDqq 

y —UJ IjJ —2X j \ (Jpq{po) 
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/-2(A-/x) 2a; \ / muaqq{oo) 

( ij^^pp I = - -2{X + ^i) -2u I I ^app(oo) I , 



then 

Dqq = (A - //)(7qq(00) Crpq(oo), 

Dpp = (A + //)(jpp(oo) + ma;^(jpq(oo), (3-21) 

1 1 

Dpq = -{rmjPaqq{oo) crpp(oo) + 2Acrpq(oo)). 

Zi TYli 

Hence, from (3.4) the fundamental constraints on a"qq(oo), a"pp(oo) and (7pq{oo) foUow: 

Dqq = (A - IJi)aqq{00) - — CTpq(oo) > 0, 

Dpp = (A + /x)(jpp(oo) + mw^(Tpq{oo) > 0, 

DqqDpp - Dlq = (A^ - ll^)(Jqq{Qo)(Jpp{oo) - U;VJg(oo) + 

+ (A - ii)muj^aqq{oo)apq{oo) - (-^ + ^) app(oo)apq(oo)- 

~{mw'^falq{Qo) - -^crlpioo) - AV2^(oo) + ^a;Vgq(oo)(7pp(oo)- 

A A^^^ 

-??7^^A(Tgq(oo)(7pq(oo) H app{oo)(Tpq{oo) > — — . (3.22) 

777' 4n 

The constraint (3.22) can be put in a more clear form: 

4(A^ +uj'^ - ii'^){aqq{oo)app{oo) - apq{oof)- 

— {muPaqq{cio) H (Tpp(oo) + 2fj,apq{oo))^ > h^X^. (3.23) 

If fj, < uj (the underdamped case), then A^ + a;^ — //^ > A^. If > a; (the overdamped 
case), then < A^ + a;^ — /i^ < A^ (A > z/) and the constraint (3.23) is more strong than 
the uncertainty inequality a"gg(oo)(Tpp(oo) — 0-^^(00) > /4. By using the fact that the 
linear positive mapping B(H) C defined by ^4 — > Tr{pA) is completely positive (hence 
2-positive), in [17] the following inequality was obtained: 



Dqqappit) + Dppaqq{t) - 2Dpqapq{t) > 



2 

From this inequality which must be valid for all values of t e (0, 00) it follows that 



Dqq(Tpp{oo) + Dppaqq{oo) - 2Dpqapq{oo) > 



n^x 



2 

Using Eq. (3.21) this inequality is equivalent with the uncertainty inequality 

aqq{oo)app{oo) - a-pg(oo) > — . 
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A restriction connecting the initial values crqq(O), crpp(O), crpq(O) with the asymptotic values 
crpp(oo), (Tqq(oo), crpq(oo) is also obtained: 



Dqqapp{0) + Dppaqq{0) - 2Dpqapq{0) > 



2 

More explicitly 

A(crgq(oo)crpp(0) + (7pp(oo)(7gq(0) - 2crpq(oo)(Jpq(0))- 

1 

-/x(cTqq(oo)cJpp(0) - (Tpp{oo)aqq{Q)) (o-pQ(oo)cTpp(O) - CTpp (oo)cTpq (0)) + 

777' 

+ma;^(cTpg (00)0-35(0) - crqq(oo)cTpq(0)) > (3.24) 
If the asymptotic state is a Gibbs state 

Pg(oo) = e jTrie fc^)^ 

then 

^^^^^^-^ " 2^ 2fcT' " coth ^,c7p5(oo) = (3.25) 

and 

A + ^ huj X — 11 % ^ ?iu> ^ , 

-Ddo = nma;coth -—-=, Dan = coth D„n = ((3.26) 

and the fundamental constraints (3.4) are satisfied only if A > // and [15]: 

If the initial state is the ground state of the harmonic oscillator, then 

, . h , . mhu) 

Then (3.24) becomes 

AfcTogfoojma; H — — ) — fiiaqqioojmuj — j > nA. 

muj muj 

For example, in the case (3.25), this implies coth.hu)/2kT > 1 which is always valid. 

Now, the explicit time dependence of aqq{t), app{t), (Tpqit) will be given for both under- 
and overdamped cases. From Eq. (3.19) it follows that in order to obtain this explicit time 
dependence it is necessary to obtain the matrix elements of Te^^T. In the overdamped 
case {ji> uj)^v^ = ji'^ — uP' we have 



Te^^T = I asi 022 023 



2\t I «11 «12 Oi3 

021 ^22 ^23 

^31 ^32 033 
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with 

011 = (a*^ + ^'^) cosh2i/t + 2(Uzv sinh 2z>'t — a;^, 

012 = (a*^ — '^^) cosh2z/t — cu^, 
ai3 = 2a;(/xcosh2zyt + i^sinh2z^t — //), 

021 = (/i^ — J^^) cosh2z>'t — a;^, 

022 = li^'^ + i''^)cosh2ut - 2/j,usmh2i't - u)'^, (3.27) 

023 = 2a;(// cosh2i/t — z/sinh2z/t — (u), 

031 = —Lv{ncosh2ut + usmh2ut — fx), 

032 = —u> {/J, cosh 2^1 — ly smb. 2^1 — /x), 

033 = -2{(jj^ cosh.2i't - jj?'). 

In the underdamped case {ji < uj),Q?' = u"^ — ji^ we have 

-2\t I bll 612 &13 

Te^*T=-^^ 621 622 623 

\ ^31 ^32 ^33 

with 

hii = (/x^ - ^2) cos 2nt - 2/xO sin 2nt - u? , 
bi2 = IfJ^ + ^ ) cos 2f]t - 

613 = 2ci;(/i cos 2nt — nsin2ni — 
621 = {fi^ + O^) cos20t - c^^ 

^>22 = (/U^ - COS 20t + 2^0. sin 20t - 0;^, (3.28) 
623 = 2u) (/X cos 2Vtt + Vt sin 2Qt — /j,), 

631 = — a;(/xcos2f2i — f2sin2f2t — /x), 

632 = — <^(a* cos 20t + O sin 2f2t — fx), 

633 = -2(a;2cos21]t - fj,^). 



4. The method of the characteristic function 

Instead of solving the master equation (3.7) directly, we first introduce the normally 
ordered quantum characteristic function x(-^) ^*jt) defined in terms of the density operator 
P by 

x(A,A*,t) =Tr[p(t)exp(Aa+)exp(-A*a)], (4.1) 

where A is a complex variable and the trace is performed over the states of system. Sub- 
stituting Eq.(4.1) into the master equation (3.7) and using the operator relations 

aexp(Aa"'") = exp(Aa"'')(a + A), 

a+ exp(-A*a) = exp(-A*a)(a+ + A*), 
Tr[p{t) exp(Aa+) exp(-A*a)(a+ + A*)] = OaX, 
Tr[p(t) exp(Aa+) exp(-A*a)a] = -dA*X 

or applying the rules: 
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+ 9 
pa -^X, 

the following partial differential equation for x is found: 

{dt + [(A - iuj)A + iiA*]dA + [(A + iuj)A* + M]5a*}x(A, A*,t) = 

= {L|A|2 + CA^ + C*A*2}x(A, A\t), (4.2) 

where 

L = X-D2,C=^{fi + Dt). 

We consider the state of the system initially to be a superposition of coherent states. The 
coherent states |q! > of the harmonic oscillator are minimum uncertainty states having 
mean coordinate < q > and mean momentum < p > given by 

/ 2fi 

< q >=< a\q\a >= \ Rea, < p >=< a\p\a >= V2hmu!lma. (4-3) 

V mu> 

Consequently, we take as the initial density operator 

p(0) = J dadpN{a,P)\a>< P\. 

The quantum characteristic function corresponding to the operator \a >< (3\ is given from 
Eq. (4.1) by 

X=< P\a> exp{XP* - X*a). (4.4) 
We look for a solution of (4.2) having the exponential form 



X 



(A,A*,t) = J dad(3N{a,P) < P\a > e^p[A{t)A + B{t)A* + f{t)A^ + r{t)A*'^ + h{t)\A\^]. 

(4.5) 

The form of the solution (4.5) is suggested from the fact that the left-hand side of Eq. 
(4.2) contains first-order derivatives with respect to the time and variables A and A* and 
is symmetric with respect to complex conjugation. The functions A{t), B{t), f{t) and h{t) 
depend only on time. Corresponding to the initial factor in Eq. (4.4), these functions have 
to satisfy the initial conditions 

^(0) = r , B{0) = a, /(O) = 0, hiO) = 0. (4.6) 
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When we introduce the function (4.5) into Eq. (4.2) and equate the coefficients for equal 
powers of A and A*, we get the following two systems of linear differential equations of 
first order with constant coefficients: 

"^^^^^ + (A - iu;)A{t) + iJiBit) = 



dt 

dB{t) 
dt 

dR{t) 



+ IJiA{t) + {X + iuj)B{t) = Q (4.7a) 



dt 

dl{t) 
dt 

dh{t) 



+ 2\R{t) + 2ujl{t) + iJih{t) = ReC 

+ 2XI{t) - 2ujR{t) = ImC (4.76) 
+ AiiR{t) + 2\h{t) = L, 



dt 

where R{t) = Ref{t),I{t) = Imf{t) with the initial conditions R{0) = /(O) = h{0) = 0. 
Subject to the initial conditions (4.6), the homogeneous system (4.7a) has the solution 
[19]: 

A{t) = u{t)(3* - v{t)a, 
B{t) = -u\t)a + v{t)l3\ (4.8) 

where 

«(0 = ^[exp(-//-t) + exp(-//+t) + — — (exp(-//+t) - exp(-//_t))], 

v{t) = (exp(-//_t) - exp(-//+t)). (4.9) 

The eigenvalues are given by 



|I± = X±^/|J,^-u^,^=^/|J,''-u^. (4.10) 



The system (4.76) has the eigenvalues — 2A, — 2(A ± \J [i^ — w^) = — 2//± and in order to 
integrate it we apply the same method as for the system (3.13) in the preceding Section. 
We obtain: 

P N iuM 

f{t) = ^ exp(-2//+t)(7 -^^)- ^ exp(-2//_t)(7 + iuj) - exp(-2At) + /(oo), 

h{t) = Mexp(-2At) + iVexp(-2//_t) + Pexp(-2//+t) + h{oo). (4.11) 
Here M, N, P, /(oo) and h{oo) are constants given by [19]: 

M = ^ililmC+^uL), 
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^ i^ReC + ujImC+^) 



27^(A + 7) 2 
and the asymptotic values connected with the diffusion coeflBcients Dgg, Dpp and Dpg are: 

2{XReC - ulmC) - L// 



R{oo) 



4(A2 - 72) 



2uXReC + 2(X'-f/')ImC-Lfiu 
'^"^^ = 4A(A2 - 72) ' ^"-^'^ 

, , , L(A2 + uj^) - 2^{XReC - ojImC) 
^^^^ ^ 2A(A2 - 72) • 

By knowing the characteristic function (4.5), (4.8)-(4.11) corresponding to the initial 
density operator which represents a superposition of coherent states, it is easy to obtain 
explicit formulae for the moments: 

< a+-^{t)a-{t) >= Tr[a+-^{t)a-{t)p{t)] = (-1)'^ g^..g^^ x(A, A*, ^)|a=a-=o- 
In the following, we take the density operator p in the coherent state representation 

p(0) = / P{a)\a >< a\(fa, 



where -P(a) is the diagonal or P Glauber distribution and d?a — dReadlma. The inte- 
gration covers the entire complex a plane. Then the characteristic function (4.5) becomes: 

X(A, A*, t)^ J d'^aP{a) exp[(ua* - va)A + {-u*a + va*)A*] exp[/A2 + /*A*2 + h\A\^]. 

Let us assume that the damped oscillator is at t = prepared in a pure coherent state, 
say \ao >, corresponding to P{a) = S{Rea — Reao)S{Ima — Imao). One has 

X^^^ (A, A*,t) = exp[{ua*o - vao)A + {-u*ao + va*o)A*] exp[fA^ + f*A*^ + h\A\'^]. 

The first moments are given by 

< a^[t) >= — — — |a=a*=o = uaQ - vao, 

< ait) >= |a^a*=o = u ao- vaQ. 



OA* 
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Then, with the notations (3. 9), using (4.3) and the transformations 



for the displacement operator q{t) and the momentum operator p{t) of the oscillator, we 
obtain the following mean values: 

with u, V given by (4.9), (4.10). There are two cases: 

a) the overdamped case: fj, > uj , = fj? — uj^ , ^ = v; then 

u{t) = exp(— At)(coshz>'t + — sinh^-t), 

v{t) = -- exp(-At) sinhi^t (4.13) 

and aq{t),ap{t) take the previous form (3.15); 

b) the underdamped case: < a;, = a;^ — 7 = ifl; then 

u{t) = exp(— At)(coshrit + ^ sinhilt), 

v{t) = exp(-At) sinhnt (4.14) 

and (Tq{t), (Tp{t) take the previous form (3.16). 
For the variances one finds: 

< a\t) ^^J^ U^A.^o = («*«o - ^a^)' + 2/*, 

< a+2(t) >= ^^^^'V =A*=o = (««o - ^«o)' + 2/, (4.15) 

< a+(t)a(t) >= — gAM* 1^^^*^° ^ ^^^o " t'Q;o)(w*Q;o - va^) - h. 

Then the relations (3.9) will give us the explicit time dependence of the variances aqq{t), 
app{t), apqit). The asymptotic values of these variances are given by the following expre- 
sions [19]: 

h 1 

CTqqioo) = (/ + /* -h + -)|t^oo, 
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crpp(oo) = -hmcvif + f* + h- -)|t^oo, 

O'pqioo) = in{f - /*)|t^oo- 

With /(oo) = R(oo) + iI{oo) and using the formulas (4.12) for i?(cxD), /(oo), /i(oo), the 
asymptotic values of the variances take the same form (3.20) as in the preceding Section, 
as expected. 

With the relations (4.15), the expectation value of the energy operator can be calcu- 
lated: 

E{t) = nu{< a+a > +^) + ^(< a+^ > - < a^ >). 
The asymptotic mean value of the energy of the open harmonic oscillator is: 

E{oo) = -^a-pp(oo) + ^mu'^aggioo) + //crpq(oo), 

or, as a function of diffusion coefficients, by using (3.20): 



5. Quasiprobability distributions for damped harmonic oscillator 

The methods of quasiprobabilities have provided technical tools of great power for 
the statistical description of microscopic systems formulated in terms of the density op- 
erator [40-42,58]. The first quasiprobability method was that introduced by Wigner [43] 
in a quantum- mechanical context. In quantum optics the P representation introduced by 
Glauber [44,45] and Sudarshan [46] provided many practical applications of quasiprobabil- 
ities. The development of quantum-mechanical master equations was combined with the 
Glauber P representation to give a Fokker-Planck equation for the laser [47,48] . One use- 
ful way to study the consequences of the master equation for the one-dimensional damped 
harmonic oscillator is to transform it into equations for the c-number quasiprobability dis- 
tributions associated with the density operator. The resulting differential equations of the 
Fokker-Planck type for the distribution functions can be solved by standard methods and 
observables directly calculated as correlations of these distribution functions. However, 
the Fokker-Planck equations do not always have positive-definite diffusion coefficients. In 
this case one can treat the problem with the generalized P distribution [41]. 

First we present a short summary of the theory of quasiprobability distributions. For 
the master equation (3.7) of the harmonic oscillator, physical observables can be obtained 
from the expectation values of polynomials of the annihilation and creation operators. The 
expectation values are determined by using the quantum density operator p. Usually one 
expands the density operator with the aid of coherent states, defined as eigenstatcs of the 
annihilation operator: a\a >— a\a >. They are given in terms of the eigenstates of the 
harmonic oscillator as 

\a\ . sr^ 1 



\a >— exp( — ) —=a'^\n > 



21 



with the normalization | < f3\a > P = exp(— |q;— In order to solve the master equation 
(3.7) we represent the density operator p by a distribution function over a c-number phase 
space. The chosen distribution function, introduced in [49] is defined as follows: 

$(ct, s) = -L y -^(A, s) exp(Q;A* - a*A)(fA, (5.1) 

with the characteristic function 

x{Ks)=Tr[pD{A,s)], 
where -D(A, s) is the displacement operator 

D{A, s) = exp(Aa+ - A*a + ^s|A|^). 
The interval of integration in Eq. (5.1) is the whole complex A plane. Because of 

S'^{a) = ~^ J exp(ctA* - Q;*A)ii^A, 

the characteristic function $(«,«) is the Fourier transform of the characteristic function. 
Since the density operator is normalized by Trp = 1, one obtains the normalization of 



/ 



In this paper we restrict ourselves to distribution functions with the parameters s = 1, and 
— 1. These distribution functions can be used to calculate expectation values of products of 
annihilation and creation operators. For that purpose one first expands the displacements 
operator in a power series of the operators a and a'^: 

D(A,.)=£5:^4=^{a+'»a"}.. (5.2) 

The braces with the index s indicate the special representations of the polynomials de- 
pending on s. For example for n = m = 1 we have the s-ordered operators: 

{a~^a}s=i — a'^a, 

{a+a}s=o = ^("^"^a + aa'^), (5.3) 
{a^a}s=-i = aa'^. 

Expectation values of the s-ordered operators can be calculated as follows: 



< {a+V}s >= Tr[p{a+VU = Tr[p(|^)-(-^)-i^(A, s)|a=o] = 
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For the last step we apply the inverse relation to Eq. (5.1): 

x(A,s) = J ^{a,s)ex.p{Aa* - A*a)(fa. 

In the following we discuss the distribution functions for s = 1,0 and —1 in more detail. 
For s = 1 we obtain the Glauber P function [44,45,50], for s = the Wigner function [43] 
and for s = — 1 the Q function [49]. For s = 1 we have 

D{A, 1) = exp(Aa+) exp(-A*a). 

Then the s ordering in Eq. (5.2) corresponds to normal ordering. Since the Glauber P 
function is the Fourier transform of the characteristic function 

Xiv(A) = Tr[pexp(Aa+) exp(-A*a)] = x(A, 1), 

it follows from Eq. (5.1) that the distribution $(a, 1) is identical to the P function. 
This function is used for an expansion of the density operator in diagonal coherent state 
projection operators [44-46,51]: 



p = J P{a)d 



a\a >< a\. 



Calculating the expectation value of normally ordered operator products we obtain the 
relation 

< a+'"a'* >= j («*)'"«"$(«, l)d^a = j {a*)"'a''P{a)d^a, 

from which we again derive P{a) — ^{a,l). Despite the formal similarity to averaging with 
a classical probability distribution, the function P{a) is not a true probability distribution. 
Because of the overcompleteness of the coherent states, the P function is not a unique, 
well-behaved positive function for all density operators. 

Cahill [52] studied the P representation for density operators which represent pure 
states and found a narrow class of states for which the P representation exists. They can 
be generated from a particular coherent state |q; > by the application of a finite number of 
creation operators. For example, for the ground state of the harmonic oscillator it is easy 
to show that xjv(A) = 1 for all A. In that case the P function becomes -P(ct) — d'^{a). The 
delta function and its derivatives are examples of a class of generalized functions known 
as tempered distributions [50] . Also Cahill [53] introduced a representation of the density 
operator of the electromagnetic field that is suitable for all density operators and that 
reduces to the coherent state P representation when the latter exists. The representation 
has no singularities. 

Sudarshan [46] offered a singular formula for the P representation in terms of an 
infinite series of derivatives of the delta function. From the mathematical point of view, 
such a series is usually not considered to be a distribution function [50,51]. 
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For s = — 1 we have 



D{A, -1) = exp(-A*a) exp(Aa+). 



The s ordering corresponds to antinormal ordering. Because the Q function is the Fourier 
transform of the characteristic function 



it follows from Eq. (5.1) that the distribution —1) is the Q function. It is given by 
the diagonal matrix elements of the density operator in terms of coherent states: 



Though for all density operators the Q function is bounded, non-negative and infinitely 
differentiable, it has the disadvantage that not every positive Q function corresponds to 
a positive semidefinite Hermitian density operator. Evaluating moments is only simple in 
the Q representation for antinormally ordered operator products. 

For s = 0, the distribution ^>(A,0) becomes the Wigner function W. The latter 
function is defined as the Fourier transform of the characteristic function 



Because this characteristic function is identical to x(A,0), we conclude that 0) is 
the Wigner function W{a). Therefore, the Wigner function can be used to calculate 
expectation values of symmetrically ordered operators: 



The symmetrically ordered operators are the arithmetic average of {m + n)\/ (mini) differ- 
ently ordered products of m factors of and n factors of a. An example for m = n = 1 
is given in Eq. (5.3). 

The Wigner function is a nonsingular, uniformly continuous function of a for all 
density operators and may in general assume negative values. It is related to the density 
operator as follows: 



Xa(A) = Tr[p exp(-A* a) exp(Aa+)] = x(A,-l), 



Q{a) = — < a\p\a > . 



TT 



Xs(A) = rr[pexp(Aa+ - A*a)] = x(A, 0). 





Also it can be obtained from the P representation: 
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By using the standard transformations [54,55], the master equation (3.7) can be trans- 
formed into a differential equation for a corresponding c-number distribution. In this Sec- 
tion we apply these methods to derive Fokker-Planck equations for the before mentioned 
distributions: the Glauber P, the Q and Wigner distributions. Using the relations 



dA 



= [(s-l)^ + a+]^(A, s) = D{A, s)[{s + l)— + a- 



= [(. + 1)^ - a]D{A, s) = D{A, s)[{s - 1)^ - a], 

one can derive the following rules for transforming the master equation (3.7) into Fokker- 
Planck equations in the Glauber P{s = 1), the Q{s = —1) and Wigner W{s = 0) repre- 
sentations: 

, s -1 d 
ap ^ (a — - — 

, . ^ s + 1 d 
pa <-> (q! — - — )<P 

Applying these operator correspondences (repeatedly, if necessary), we find the fol- 
lowing Fokker-Planck equations for the distributions $(q;,s): 

Here, $(«,«) is P{s — l),Q{s — —1) or W{s = 0). While the drift coefficients are the 
same for the three distributions, the diffusion coefficients are different: 

da = -(A -I- iuj)a + fia*, Daa = -Di + s/i, D^a* = D2- sA. 

The Fokker-Planck equation (5.4) can also be written in terms of real coordinates xi 
and X2 defined by 



' muj 1 
a = xi + ix2 = \l -rrr < q > +i , < P>, 



a* = xi — 1X2 (5-5) 



as follows: 
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with the new drift and difusion coefficients given by 

di = -(A - IJ,)xi + UJX2, d2 = -LOXi - (A + /x)x2, 



We note that the diffusion matrix 



D = 



1^11 1^12 
Dl2 D22 



for the P ddistribution (s = 1) needs not to be positive definite. 

Since the drift coefficients are hnear in the variables xi and X2{i = 1, 2): 



2 

di — ^ ] ■^ij'^ji -^ij — 



ddj 
dxj 



with 



1 (5.7) 



and the diffusion coefficients are constant with respect to xi and X2, Eq. (5.6) describes 
an Ornstein-Uhlenbeck process [56,57]. 

The solution of the Fokker-Planck equation (5.6) can immediately be written down 
provided that the diffusion matrix D is positive definite. However, the diffusion matrix in 
the Glauber P representation is not, in general, positive definite. For example, if 

DfiDfa - {Di2f < 0, 

the P distribution does not exist as a well-behaved function. In this situation, the so-called 
generalized P distributions can be taken that are well-behaved, normal ordering functions 
[41]. The Q and W distributions always exist; they are Gaussian functions if they are 
initially of Gaussian type. 

From Eq. (5.6) one can directly derive the equations of motion for the expectation 
values of the variables xi and X2{i = 1,2): 

d < Xi > y-^ . , ^ , 

^ = -E^^^- <^j>- (5-8) 

By using Eqs. (3.5), (5.5) and (5.8) we obtain the equations of motion for the expec- 
tation values aq{t), ap{t) of coordinate and momentum of the harmonic oscillator which 
are identical with those derived in the preceding two Sections by using the Heisenberg 
representation and the method of characteristic function, respectively (see Eqs. (3.12)). 
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The variances of the variables xi and X2 are defined by the expectation values 



They are connected with the variances and covariance of the coordinate q and momentum 
P by 

o'qq = {2h/mu!)aii, app = 2hmuja22, 

^pq ^{pq + qp) > - < p >< q >= 2hai2. 

They can be calculated with the help of the variances of the quasiprobability distributions 
(*,J = 1,2): 

The following relations exist between the various variances: 

(^ii - (^ii + ^ - (^ii - ^ - (^ii - 2, 

^ — „P — „Q — „w 
(^12 — cr 12 — cr -^2 — ^12- 

The variances cr^*^ fulfill the following equations of motion: 

-if- = -j:iAuaif+-ii'^) + Dff. (5.9) 
They can be written explicitly in the form: 



^ = -2A,,aif-2A,2aii + D[\\ 



dt 

= —2^210-12^ — 2^220-22^ + -D22^ 

J (s) 

— -{All + A22)(Ji2 -A2i(Jii - Ai2a22 + J^12 i 

where the matrix elements Aij are defined in Eq. (5.7). These relations are sufficient to 
prove that the equations of motion for the variances an and (722 and the covariance a 12 
are the same irrespective of the choice of the representation as expected. The correspond- 
ing equations of motion of the variances and covariance of the coordinate and momentum 
coincide with those obtained in the preceding two Sections by using the Heisenberg repre- 
sentation and the method of characteristic function, respectively (see Eqs. (3.13)). 
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In order that the system approaches a steady state, the condition X> i/ must be met. 
Thus the steady-state solutions are 

<^{xi,X2,s) = - — , exp[-- J] (cT-^)y(oo)a:ia:j], (5.10) 

2TTy/ det{a [oo)) ^ . 2 

where the stationary covariance matrix 

.(co)=.WM=(ii;(-' 

V (^12 (00) ay {00) J 
can be determined from the algebraic equation (see Eq. (5.9)): 

1=1 

With the matrix elements Aij given by (5.7), we obtain 

(s). , _ (2A(A + u;^)d[]' + uj^Pji + 2a;(A + //)£>g 
"^11 (°°) — 



4A(A2 + a;2 _ ^2) 

(7^(00) 



_ ^'^ff + (2A(A - //) + a;^)i^g - 2a;(A - f,)D[i 



4A(A2 + a;2 - //2) 

(,) -a;(A + ii)d[^ + c^(A - ^)Dg + 2(A^ - /.^)i^g 
= 4A(A2 + c.2-^2) • 

The explicit matrix elements cr^^ for the three representations P, Q and W can be obtained 
by inserting the corresponding diffusion coefficients. The distribution functions (5.10) can 
be used to calculate the expectation values of the coordinate and momentum and the 
variances by direct integration. The following relations are noticed [37] : 

<^^{<^) = l^ioo) + = 1, 2. 

The uncertainty principle crii(722 > 1/16 gives rise to the conditions 4cr^(T^ > + cr^ 
and (y\\(y^2 > 1/16 for the Q and W distributions, respectively. 



6. Density matrix of the damped harmonic oscillator 

In this section we explore the general results that follow from the master equation 
of the one-dimensional damped harmonic oscillator. Namely, we discuss the physically 
relevant solutions of the master equation, by using the method of the generating function. 
In particular, we provide extended solutions (including both diagonal and off-diagonal 
matrix elements) for different initial conditions. 
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The method used in this section follows closely the procedure of Jang [39]. Let us 
first rewrite the master equation (3.7) for the density matrix by means of the number 
representation. SpecificalUy, we take the matrix elements of each term between different 
number states denoted by |n >, and using a~^\n >— \/n+ l\n+l > and a\n >— y/n\n — l > 
we get 

= -i(jj{m - n)pmn + >^Pmn - (m + H + l)D2Pmn + 

Pm— l,n+l~l~ 

1 1 

+ -^J{n+l){n + 2){D^ - fJ,)pm,n+2 + -y/{m+ l){m + 2){Dl - fJ,)pm+2,n- 

-^/rnTl)nDlprn+l,n-l + ^y^in^^T)n{Dl + At)Pm,n-2+ 

+ V^(rrr+lJ(Tr+l)(i:'2 + A)p^+i,n+i + y/mn{D2 - X)pm-i,n-i- (6.1) 
Here, we have used the abbreviated notation 

Pmn =< m\p{i)\n > . 

This master equation is complicated in form and in indices involved. It comprises not only 
the density matrix in symmetrical forms, such as Pm±i,n±i, but also those matrix elements 
in asymmetrical forms like Pm±2,n, Pm,n±2 and Pm^i,n±i- In order to solve Eq. (6.1) we 
use the method of a generating function which allows us to eliminate the variety of indices 
m and n implicated in the equation. When we define the double-fold generating function 

by 

G{x, y,t) = J2 -7^x"^y^pmn{t), (6.2) 
^ Vm!n! 

the density matrix can be evaluated from the inverse relation of Eq. (6.2): 

Pmnit) = -yl={^r{^rGix,y,t)U=y=o, (6.3) 

vmini ox ay 

provided that the generating function is calculated beforehand. When we multiply both 
sides of Eq. (6.1) by x^y"^ j \J m!?i! and sum over the result, we get a linear second order 
partial differential equation for G(x, t), namely 

— G(a;, 2/, t) = {[-{iu + D2)x - Dly]— + [-Dix + {iu - ^2)2/]^ + 

+ [^{D, + p)x^ + ^-{Dl+p)y^ + {D2 - X){xy - l)]}G{x,y,t). (6.4) 
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A special solution of Eq. (6.4) can be taken as 

Gix,y,t) = i exp{xy - [B{x - Cf + L>(y - Ef + F{x - C){y - E)]/H}, (6.5) 

where A, B, C, D, E, F and H are unknown functions of time which are to be determined. 
When we first substitute the expression (6.5) for G{x, y, t) into Eq. (6.4) and equate the 
coefficients of equal powers of y and xy on both sides of the equation, we get the following 
differential equations for the functions B, D, F and H: 

-J^ = -i^l - H)^ -{D,-n)^- {D, + A)^ + 2A, (6.6) 

|(|) = 2(A - ia;)| - /.^ - \{D, - /^)^ - 2{D^ + A)^ - 2{Dl - ^)|^, (6.7) 

-(-) = 2(A + za;)- - ^- - -{D* - - 2{D, + A)-^ - 2{D, - ^)-^, (6.8) 

fj w z? p2 f} D H F{ J~) E 

^(-) = 2A--(D.+A)^-2,(^ + ^)-4(0.+A)^-2(D;-,)^-2(D,-,)^. 

(6.9) 

In addition to these equations, we get for the functions C and E 

2B^ + F^ = (-2(A - iuj)B + i^F)C + (2/xS - (A + iuj)F)E, (6.10) 

2D^ + F^ = (-2(A + iuj)D + i^F)E + (2yuL> - (A - iLo)F)C. (6.11) 

The equations (6.10) and (6.11) can be reformulated in order to eliminate the functions 
B, D and F, provided BD - F^/A ^ 0. We obtain 

dC 

= -(X-iu;)C + nE, (6.12) 

dE 

= -(^X + iu;)E + nC. (6.13) 

(JjL 

The functions A, B, D, F and H are connected by the auxiliary condition that Trp is 

independent of time. The trace of p can be evaluated by summing the diagonal matrix 
elements pnn given in Eq. (6.3) or directly by using the integral expression 

Trp — Pnn = 77—^ / exp(-/ci/c2) exp(zA;ia; + ik2y)G{x, y, t)dkidk2dxdy . 
^0 ^27r)^ J 

We obtain with the generating function (6.5) 

Tr{p) = {^{^ - BD))-^'\ (6.14) 
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This quantity is time-independent which can be verified by constructing an equation sat- 
isfied by the quantity {F'^/A - BD)/H^. Combining Eqs. (6.7)-(6.9) we get 

We see immediately that the first factor on the right-hand side of this equation is identical 
with the right-hand side of Eq. (6.6). Accordingly, we find 

Since the scaling function H is arbitrary, we simplify the following equations by the choice 

p2 

— -BD = -H. (6.15) 

Setting Trp = 1, we obtain from Eqs. (6.14) and (6.15) the normalization constant = 
—if/4. As a consequence, we can simplify Eqs. (6.7)-(6.9) by eliminating the function H 
from these equations. The resulting three equations are 

^ = -2(A + iu;)B -i^F + 2{Di - //), 

^ = -2{X-iuj)D - nF + 2{D*- n), (6.16) 
dF 

— = -2n{B + D)- 2XF - A{D2 + A). 

These equations imply that the function D is complex conjugate to B, provided that the 
function F is real. 

In order to integrate the equations for the time-dependent functions B,C,D,E and 
F we start with Eqs. (6.12) and (6.13). These equations imply that the function E is 
complex conjugate to the function C. By solving the coupled equations we find: 

C{t) = E*{t) = u{t)C{0) - v{t)C*{0), (6.17) 

where u{t) and v{t) are given by (4.13) and (4.14) for the two considered cases: overdamped 
and underdamped, respectively. For integrating the system (6.16) we proceed in the same 
way as for integrating the system (3.13). With the assumption that F is real and 

D{t) = B*{t) = R{t)+il{t), 

we obtain explicitly: 

R{t) = ^{e-^^+' + e-^^-')R + ^(6-2/^+* - e-^^-')(^I+ ^F) + R{oo), 
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nt) = -e--(^/ + ^F)+ 



7^ 7 



where we used the notations: 



//± = A±7,7= 



R = R{0) - i?(oo), / = /(O) - /(oo), F = F(0) - F(oo). 

We can also obtain the connection between the asymptotic values of B{t), D{t), F{t) and 
the coefficients Di, D2, A* 

A(i?eL>i - //) + a;/mL»i + //(L>2 + A) 
i?(oo) = ReDioo) = — — ^ -, 



/(oo) = /m-D(oo) 
F{oo) = -2 



A^ -7^ 

ujX{ReDi - + iii'^ - X^)ImDi + a;/i(i:>2 + A) 
A(A2 - 72) ^ 

//[A(i?ei:>i - //) + ujImDi] + (A^ + a;2)(i:)2 + A) 



A(A2 - 72) 



When all explicit expressions for B, C, D, E, F and i/ are introduced into Eq. (6.5), we 
obtain an analytical form of the generating function G{x, y, t) which allows us to evaluate 
the density matrix. 

If the constants involved in the generating function satisfy the relations 

C(0) = 0, R{0) = i?(oo), 7(0) = /(oo), F(0) = F(oo), 

we obtain the stationary solution 

C{t) = E{t) = 0,R{t) = R{0),I{t) = I{0),F{t) = F(0), 

so that 

D{t) = B*{t) = i?(0) + z/(0), 
H{t) = -4A^{t) = R^{0) + /^(O) - F2(0)/4. 
Then the stationary solution of Eq. (6.4) is 

Gix, y,t) = j exp{(l - ^)xy - {Bx^ + B*y^)/H}. (6.18) 
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In addition, for a thermal bath [17] with 



n nmw 

the stationary generating function is simply given by 

The same generating function can be found for large times, if the asymptotic state is a 
Gibbs state with = 0. In this case we obtain with Eq. (3.26) and = 

D2 = A coth 



2kT 

and ^ 

G{x, y) = (1 - exp(-^)) exp{exp{-—)xy). 

The density matrix can be calculated with Eq. (6.3) and yields the Bose-Einstein distri- 
bution 

p = (1 - exp(-^)) exp{--j-^)Snm- 

A formula for the density matrix can be written down by applying the relation (6.3) 
to the generating function (6.5). We get 



< m\p{t)\n >= ^— J- exp[-{BC^ + DE^ + FCE)/H]x 

^ M _ iL\n3 (B_\ni (D_\n2 (2^^ 4- PE \m-2n^ /o DE , FC_\n-2n2—nz 

\ V H ' ^ H ' ^ H ' ^ H H ' ^ H H ' 1 Q"\ 

^ „ ni!n2!n3!(m — 2ni — n3)!(n — 2n2 — ns)! 

ni,n2,n3=0 

In the case that the functions C{t) and E{t) vanish, the generating function has the form 
of Eq. (6.18). Then the elements of the density matrix with an odd sum m + n are zero: 
Pmn = 0form + n = 2/c-|-l with k = 0, 1, 2, ... The lowest non- vanishing elements are 
given with pmn = Pnm as 

_ 1 _ V2B _ 1 i", 

P00-^,P20--^,Pll-^(l-^), 

2BB* 1, F 2 F \fQB VGB"^ 

It is also possible to choose the constants in such a way that the functions B and D vanish 
at time t = and F(0) = H{0). Then the density matrix (6.19) becomes at t = 0{E = C*): 

< m|p(0)|n >= ^l=(C*(0))-(C(0))"exp(-|C(0)|2). (6.20) 
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This is the initial Glauber packet. The diagonal matrix elements of Eq. (6.20) represent a 
Poisson distribution used also in the study of multi-phonon excitations in nuclear physics. 
In the particular case when we assume 

Di= ii = 0,D2 = X, 

B(0) = D{0) = 0, F(0) = H(0) = -4, 

the differential equations (6.16) yield B{t) = D{t) = and F{t) = H{t) = -4. Then the 
density matrix subject to the initial Glauber packet is (see also [39] ) 

< m\p{t)\n >^ -yl={C*{t)r{C{t)reM-\C{t)n 

where C{t) is given by Eq. (6.17). 

4. Conclusions 

The Lindblad theory provides a selfconsistent treatment of damping as a possible 
extension of quantum mechanics to open systems. In the present paper first we studied 
the damped quantum oscillator by using the Schrodinger and Heisenberg representations. 
According to this theory we have calculated the damping of the expectation values of 
coordinate and momentum and the variances as functions of time. The resulting time 
dependence of the expectation values yields an exponential damping. Second we have 
also shown how the quasiprobability distributions can be used to solve the problem of 
dissipation for the harmonic oscillator. From the master equation of the damped quantum 
oscillator we have derived the corresponding Fokker-Planck equations in the Glauber P, the 
antinormal ordering Q and the Wigner W representations and have made a comparative 
study of these quasiprobability distributions. We have proven that the variances found 
from the Fokkcr-Planck equations in these representations are the same. We have solved 
these equations in the steady state and showed that the Glauber P function (when it 
exists), the Q and the Wigner W functions are two-dimensional Gaussians with different 
widths. Finally, we have calculated the time evolution of the density matrix. For this 
purpose we applied the method of the generating function of the density matrix. In this 
case the density matrix can be obtained by taking partial derivatives of the generating 
function. The generating function depends on a set of time-dependent coefficients which 
can be calculated as solutions of linear differential equations of first order. Depending 
on the initial conditions for these coefficients, the density matrix evolves differently in 
time. For a thermal bath, when the asymptotic state is a Gibbs state, a Bose-Einstein 
distribution results as density matrix. Also for the case that the initial density matrix is 
chosen as a Glauber packet, a simple analytical expression for the density matrix has been 
derived. The density matrix can be used in various physical applications where a Bosonic 
degree of freedom moving in a harmonic oscillator potential is damped. For example, one 
needs to determine nondiagonal transition elements of the density matrix, if an oscillator 
is perturbed by a weak electromagnetic field in addition to its coupling to a heat bath. 
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The density matrix can also be derived from the solution of the Fokker-Planck equation 
for the coherent state representation. 
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